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We present a set of coupled continuum equations with a specific coupfing between mobile grains 
p and clusters h on the surface of a sandpile. The equations are analysed self-consistently; we 
demonstrate that Edwards' infrared divergence is responsible for the unexpected critical exponents 
we find, which are verified by simulations. 
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A certain class of coupled Langevin equations has been useful of late in describing the dynamics of sandpile 

^ I ^ ' surfaces [Q. The fluctuations on the sandpile surface are described by the local height h{x, t) and the density p{x, t) of 

, flowing grains - the flow being initiated by the formation of bumps. The coupling between the two variables is described 

' in terms of a transfer term which converts static grains to mobile grains and vice-versa. In this communication we 

1 ; [ point out that a singularity discovered by Edwards Q three decades ago in the context of fluid turbulence is present 
d • in such models due to the specific choice of the transfer term. This singularity largely controls the dynamics and 

' produces unexpected exponents. Our contention is supported by numerical evidence. 

I , We begin with a model introduced by Mehta, Luck and Needs (MLN) Q where the coupled Langevin equations 

' are 

O ; — ^ DhV^h~T{h,p) + 7^hix,t) 

^^DpV'p + T{h,p) + 7jp{x,t) (1) 
T{h,p) = -pp{Vh) 

where the terms rifi{x,t) and rjp{x,t) represent Gaussian white noise. A variant of the above is the model due to 
Bouchaud et al. (BCRE) |] where 
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■ r = ~vVh - pp{Vh) 

o\'. 

' and the noise is present only in the equation of motion for h. All the results which we will find here for the MLN 
I system also holds for BCRE which we have checked both analytically and numerically. 
Q ■ A simple physical picture of the coupling or 'transfer' term T{h, p) between h and p is the following: flowing grains 
are added to regions of the interface which are at less than the critical slope, and vice versa, provided that the local 
density of flowing grains is always non-zero and in proportion to the local density. This is the simplest possible form 
of exchange between the species that appears intuitively reasonable during the process of avalanching and we analyse 
Q , in what follows the profiles of both species consequent on this form. 
IL/ ■ Before this, we review some well-known facts about interfacial roughening Three critical exponents, a, /3, and 
z, characterise the spatial and temporal scaling behaviour of a rough interface. They are conveniently defined by 
considering the (connected) two-point correlation function of the heights 

H ■ 

; S{x-x',t-t') = {h{x,t)h{x',t')) - {h{x,t)){h{x'X)). 
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We have 



S'(a;,0) - (|a;|->oo) and ^(a; = 0, i) - (|i| -> c5o), 

and z = a/ (3. 

In recent years a self-consistent mode coupling analysis used hitherto in dynamic critical phenomena has been 
used to look at in particular the Kardar-Parisi-Zhang (KPZ) equation |^ and we extend its use to the case of the 
coupled equations presented here. 

In this method we set up equations (to one-loop order) for the correlation functions and self-energies in terms of the 
full Green's functions, correlation functions and vertices using assumed scaling forms for each. The critical exponents 
a and (3 defined above are obtained from the self-consistent solutions of these equations after setting Dh and Dp to 
unity. 

The analysis of these functions will be in terms of a weak scaling hypothesis which states 



along with a similar scaling relation for Gp{k, uj). A strong scaling would imply the existence of a single time scale i.e. 
Zh = Zp. As we show below, this cannot be the case here. The absence of strong scaling implies that the roughness 
exponents ah and ap may become functions of k. 

We consider the full Green's function Gh{k,tjj), which is given via the well-known Dyson equation, 

Gj^\k,u) = Gl'\k,u:) + ^n{k,u) 
Here, the zeroth order Green's function is 

Gl{k,u) = (-*C. + fc2)-l 

In the scaling limit fc^ can be dropped in comparison with Yih{k,ijj). To one-loop order, the self-energy Y^h{k) is 
given by 

T,hik,uj) = ii^{2tt)-^ J dq J dflGhik ~ q,uj - n)Spiq,n) Hk - q) (2) 

We note that due to the presence of the term Sp{q, fl), the integral is dominated by the singularity in the integrand at 
q — > 0. This 'infrared divergence' which relates to the divergence of the internal momenta q, is very different from the 
usual divergences encountered in critical phenomena where the latter occur for small wave numbers and are associated 
with long wavelength instabilities in the external momenta. In this case due to the infrared divergence in the above 
equation in the internal momenta q, the integral diverges for any value of the external momenta k, so long as ap > 0. 
This is the divergence found by Edwards for the Navier-Stokes equation We thus need either to evaluate the 
integral with a lower cut-off fco or to introduce a suitable regulator. We follow the first of these procedures for the 
above equation. 

We then proceed to evaluate the self-energy at zero external frequency, i.e Y,h{k,uj = 0) from Eq.(|2|). The integral 
in Eq.(H) becomes in the limit of zero external frequencies 

,21,2 



We have to evaluate the integral by cutting off the momentum integration at fco ^ 1 , we follow the first of the 
procedures given above to handle the infrared divergence. This gives, after some simplification. 

From the above equation with the scaling relation S^(fc) ~ fc^'* we find, on equating powers of fc, 

Zh ^ I 

We note here that the presence of the term pV/i could in principle cause the vertex to renormalize, leading to a 
correction to Zh. We have checked that this correction vanishes as fc ^ in the lowest order. 
The structure factor at one-loop level is given by 
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Sh{k,uj) 



C^2 + |I]^(fc,C^)|2 

which on integrating with respect to to gives Sh{k, t ~ 0) as 



(3) 



S^ik,t = 0)^js,ik,u)^ = ^ + ^ (4) 

Recognising that the scahng form of Sh{k,t = 0) ~ fc~^~2ah^ -^^^e notice that ah cannot in general be determined from 
Eq.(^. This is because the second term on the right-hand-side of Eq.(Q) dominates at smaU momenta k provided 
ah > 0. 

We turn now to the critical exponents in p. The single loop self-energy Y,p{k,Lu) is given by 

Ep(fc, LU = 0) = -fi^{2TT)-^ j dq j dnCpik - q, -n)Sh{q, fl)q^ (5) 

This gives, on performing the integral over internal frequency f2, 

^p{k,uj^ 0) = V / ^ i (6) 

We see from the above that Sp(fc, 0), the relaxation rate for p fluctuations, is negative and finite as fc ^ 0, and we need 
to add a positive constant, Sq, to the self-energy (Eg > |Sp(A: — > 0)|) for regulatory purposes. This divergence in the 
relaxation rate, needing regulation, is reflected in the divergence we have encountered in our numerical investigations 
below; we have there followed an analogous procedure by introducing a numerical regulator which replaces divergent 
values of the transfer term by suitably defined cutoffs j2|. The resulting constancy of Ep implies Zp « for the 
regulated equations and will be used in the following. 
The correlation function Sp{k,uj) is given by 



Sp{k,Lu) ^ {lu^ + k'''')-\2Tr)-' J dq J dn{k ~ qyShik - q,uj - n)Sp{q,n) (7) 
which on itegration over uj gives 

Sp{k,t^ 0) - fc-(l+2"p) ^ fcl-2a,+z, 1 ^ 

Finally using Zp « we have 

oip = ah + ^ — I for large k (8) 
ap ^ ah — 1 for small k (9) 

Given our numerical result of ah = 0.5, the above predicts a negative ap, at small k. This is consistent with, and 
validates our assumption of, a cutoff fco which arises naturally as the wavevector separating the region of ctp < and 
ap > 0. 

The coupled equations have been numerically integrated by using the method of finite differences. Our grids 
in time and space were kept as fine-grained as computational constraints allowed. This is in order to avoid the 
instabilities associated with the discretisation of nonlinear continuum equations. Convergence has been checked by 
keeping At small enough such that the quantities under investigation are independent of further discretisation. In all 
our calculations, we chose Dh = Dp = u = 1. 

On discretising the equations Eqs.(|^) we found once again the divergences that were previously observed in 
These divergences are in our view a direct representation of the negativity of Sp. We follow here a parallel course in 
regulating these via an explicit regulator. In earlier work ||^, a regulator was introduced which replaced the function 
ppVh by the following: 

T = +1 for ppiVh) > 1 

= pp{Vh) for - 1 < pp{Vh) < 1 
= -1 for pp{Vh) < -1 
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In addition in this paper, we have introduced noise reduction to the regulated equations which has led to a more 
accurate evaluation of all our critical exponents. 

Our results for the single Fourier transforms for the (hh) correlation function are 

(i) The Fourier transform Sh{k,t ^ 0) (Fig.|l|) is consistent with a spatial roughening exponent ^ 0.5 ± 0.02 via 
our observation of Sh{k,t^ 0) - fc-2-00±o.04 

(ii) The Fourier transform Sh{x = 0,cj) (Fig.^) is consistent with a temporal roughening exponent /3 ^ 0.48 ± 0.01 
via our observation of Sh{x = 0, ~ (j;-i-96±o.02 jjgj^gg ^ \ o.02 consistent with our prediction. 

The full structure factor Sh{k,uj) has been calculated at three different k points and Fig.^ displays a fit of our results 
to an appropriately scaled form of Eq.(||). The spatial structure factor Sh{k,u> = 0) shows a power-law behaviour 
(Fig.^) given by Sh{k, cj = 0) ~ ^-3.30±.(K qualitative accord with our expression, which predicts an exponent of —3. 
The temporal structure factor Sh{k = 0,uj) shows a power-law behaviour (Fig.^j) given by Sh{k — 0,lu) ^ ll|-1-97±-03 
which is in agreement with our prediction of cj^^ from the expression of the structure factor (Eq.(^). 

Given our values of ah — 0.5 and Zh — 1, Eqs.(^ and predict a crossover in ap from 0.0 at large k to -0.5 
as fc 0. We observe that the single Fourier transform Sp{k,t = 0) (Fig.||) shows a crossover behaviour from 
Sp{k,t = 0) ^-2.oo±o.08 large wavevectors to Sp{k,t = 0) ~ constant as fc ^ 0. This indicates a crossover from 
ttp = 0.5 for large k to -0.5 as fc ^ 0, which shows the same trend as the prediction above. Note however that the 
simulations also manifest in addition to the above the normal diffusive behaviour represented by ap — 0.5 at large 
wavevectors. This anomalous smoothing behaviour in p is the direct consequence of the infrared divergence in Eq.(0) 
discussed in the preceding. 

We have analysed in the above via perturbative methods a model of sandpile dynamics which was presented without 
analysis in earlier work The good agreement between our theoretical predictions and numerical simulations 
confirms our contention that Edwards' infrared divergence |^ plays a crucial role in producing unexpected critical 
exponents in our model. 
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FIG. 1. The single Fourier transform Sh{k,t — 0). with a power-law fit to an exponent of 2.0 ± 0.04. 

FIG. 2. The single Fourier transform Sh{x — 0,lj). with a power-law fit to an exponent of 1.96 ± 0.02. 

FIG. 3. The double Fourier transformrs Sh{ki,Lo) evaluated at three different ki = 0.1, 0.2, 0.3 fitted to Eq.(^. 

FIG. 4. The double Fourier transform Sh(k,uj = 0) with a power-law fit to an exponent of 3.3 ± 0.05. 

FIG. 5. The double Fourier transform Sh{k — 0, lo) with a power-law fit to an exponent of 1.97 ± 0.03. 

FIG. 6. The single Fourier transform Sp{k,t = 0) with a power-law fit to an exponent of 2.0 ± 0.08 at large k crossing over 
to an exponent of zero at small k. 
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